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1.  INTRODUCTION 


The  dissemination  of  solid  payload  material  as  a  gaseous  cloud  from  a  container  can 
be  accomplished  using  a  pressurized  system  (e.g.  compressed  CO2  cartridge),  a  mechan¬ 
ical  system  (e.g.  plunger)  or  a  pyrotechnic  system  (e.g.  central  burster  or  hol/cohi  gas 
flow).  In  the  latter  technique,  dissemination  of  the  payload  is  realized  after  ablation  and 
vaporization  of  the  material  from  exposure  to  the  thermal  effects  of  pyrotechnic  combustion. 
Experimental  programs  at  the  U.S.  Army  Chemical  Research  Development  and  Engineer¬ 
ing  Center  (CRDEC)  are  being  used  to  test  the  efficiency  of  the  thermal/ablation  metliod. 
The  goals  are  high  rate  of  dissemination  and  low  temperature/high  payload  concentration 
in  the  dissemination  cloud.  These  tests  can  involve  highly  instrumented  full-scale  and  small- 
scale  grenade  models  ignited  within  large  test  chambers.  Numerical  simulation  in  support  of 
such  tests  was  conducted  at  the  U.S.  Array  Research  Laboratory  (ARL).  These  simulations 
are  required  to  aid  in  data  reduction,  contribute  to  the  physical  understanding  of  thermal 
dissemination,  and  conduct  parametric  studie.  that  may  guide  future  tesls. 

The  internal  design  of  a  grenade  used  for  the  thermal  dissemination  of  solid  payload  into 
the  atmosphere  can  consist  of  two  concentric  cylinders;  a  pyrotechnic  device  in  the  outer 
annulus  and  payload  material  bonded  to  the  wall  of  the  inner  cylinder.  The  two  chambers 
are  connected.  Combustion  of  the  pyrotechnic  produces  a  high  pressure  within  the  grenade. 
A  pressure  difference  between  the  atmosphere  and  inside  the  grenade  induces  a  thru-flow 
that  thermally  erodes  and  vaporizes  the  material  in  the  inner  chamber.  The  material  in 
gaseous  form  is  entrained  in  this  flow  and  expelled  from  the  grenade.  Figure  1  shows  the 
internal  schematic  of  a  generic  grenade  configuration. 

Numerical  simulation  of  the  fluid  dynamics  associated  with  payload  dissemination  from 
a  grenade  is  accomplished  in  the  present  study  using  the  Navier-Stokes  equations  along  with 
equations  that  govern  chemical  species  ablation  and  diffusion.  Chemical  reactions  bet  ween 
payload  and  pyrotechnic  gases  may  be  but  are  are  not  necessarily  involved.  An  implicit  fiiiitc- 
difference  scheme  based  on  successive-over-relaxation  is  used  to  solve  these  equations  for  Uic 
physical  domain  of  interest.  This  domain  resides  within  the  inner  chamber  of  the  grenade; 
thus,  the  flow  within  the  grenade  is  modeled  excluding  combustion  of  the  pyrotechnic.  Using 
chamber  dimensions  and  payload  chemical  properties  this  simulation  yields  velocity,  pre.ssure, 
temperature,  density  and  chemical  compostion  of  the  gas  in  the  inner  chamber  and  exiting 
the  grenade. 
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2.  NUMERICAL  SIMULATION 


The  RAMCOMB  (RAMjet  COMBustion)  computer  code  was  originally  developed  for  the 
numerical  simulation  of  combustion  in  a  tubular  solid-fuel  ramjet  (SFRJ)  projectile  (Nusca 
et.al.  1988,  Nusca  1989).  Solid  fuel  regression  rate  and  projectile  thrust  predictions  compared 
favorably  with  in-flight  and  ground  test  data.  For  the  SFRJ  application  the  RAMCOMB 
code  simulated  a  mass- controlled  (stoichiometric)  reaction  of  non-premixed  solid  fuel  and 
oxygen  using  classical  diffusion  flame  techniques.  The  code  has  also  been  used  to  simulate 
finite-rate  premixed  gaseous  fuel  combustion  in  the  ram  accelerator  projectile  launch  system 
(Nusca  1991)  with  reaction  rates  formulated  in  terms  of  temperature  and  chemical  species 
mass  fraction.  Application  of  the  code  to  payload  dissemination  simulation  for  grenades 
involves  chemical  species  ablation  and  diffusion  without  chemical  reactions.  The  governing 
equations,  boundary  and  initial  conditions  as  well  as  the  solution  method  are  outlined  below. 


2.1  Governing  Equations.  Since  the  grenade  payload  chamber  geometry  is  axisym- 
metric  (Fig.  1)  the  governing  equations  can  be  written  in  cylindrical  coordinates.  The 
velocity  components  in  this  system  are  and  w  for  the  radial  (r),  azimuthal  (0),  and  axial 
(2)  directions,  respectively.  Axisymmetric  flow  is  aissumed  thus,  all  ^-derivatives  are  ignored; 
however,  the  azimuthal  velocity  component  and  the  azimuthal  momentum  equation  are  re¬ 
tained  (for  future  consideration  of  chamber  rotation).  Since  steady  flows  are  considered,  time 
derivatives  {djdt)  are  ignored.  The  conservation  equations  for  global  ma.ss,  momentum  (ra¬ 
dial,  axial,  azimuthal)  and  species  mass  conservation  are  given  by  (Nusca  1991,  Schlichting 
1979), 


r  dr  az 

(1) 

V  •  [puV  -  fr]  -  -Tee]  +  ~  =  0 

(2) 

V  •  [pwV  -  t"]  -I-  ^  =  0 

(3) 

-V  •  [r{pvV  -  r^)]  =  0 
r 

(4) 

d  ,  d 

—  [rpuTTij  -f  rjj^r)  +  ^  [rpwmj  + 

=  0  (5) 

Energy  conservation  for  a  compressible  flow  is  expressed  by  the  First  Law  of  Thermody¬ 
namics.  The  steady  form  of  the  First  Law  states  that  the  net  rate  of  stagnation  enthalpy 
(A  =  /i  -f  V^/2)  inflow  for  a  control  volume  is  equal  to  the  sum  of  the  shear  work  done 
by  the  contents  of  the  control  volume  on  the  surroundings  (f)  and  the  heat  transfer  to  the 
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surroundings  (Schlichting  1979). 


V- 


pVh  +  .4  +  X!]  +  Jk~  (uTt  +  vTf  + 


=  0 


(6) 


where  J  is  a  flux  term  for  mass  (Jj  =  heat  ( =  {p^lVt)Cp^VT  = 

TfcCpjVr),  and  turbulence  kinetic  energy,  k,  (J*  =  (/ieir/Pr)VA:  =  r*Vjt).  Where  F  repre¬ 
sents  the  diffusion  coefficient.  The  mass  fraction  and  molar  specific  enthalpy  for  species  j 
are  rrij  and  hj,  respectively.  The  density,  pressure,  velocity  vector  and  velocity  magnitude 
are  given  by  p,p,  VyV  ~  i/u*  +  v^  +  u?^,  respectively. 


In  Equations  2-4,6  the  shear  stress  (f)  includes  the  Reynolds  stress  with  an  effective  fluid 
viscosity  expressed  as  the  sum  of  the  molecular  and  turbulent  viscosities,  peff  =  p  +  Pt-  The 
flow  Reynolds  number,  Re,  represents  the  ratio  of  mass  flux  {pVL)  to  fluid  viscosity,  peir- 
Molecular  viscosity  (p)  is  defined  using  Sutherland’s  expression  (Ames  1958), 

'T’1.5 

p  =  2.270x10-*  - . (7) 

^  T+ 198.6  '  ' 

Turbulent  viscosity  (pt)  is  described  in  the  next  section  of  this  report. 


The  calorically  perfect  gas  assumption  csm  be  used  to  determine  the  specific  heat  of  each 
species,  Cp^,  when  the  temperature  dependence  of  the  species  is  not  well  determined.  The 
specific  heat  can  also  be  formulated  using  an  explicit  temperature  dependence  obtained  from 
tabulated  data  (Stull  and  Prophet  1971). 


Cp,/»  =  Ax  +  A{r  +  AaT*  -h  A^T^  +  A,T* 


(8) 


where  the  mixture  specific  heat  is  given  by, 

N 

<^  =  X]mjCp, 


(9) 


Mixture  temperature  (T)  is  obtained  from  the  conservation  of  energy  (Eq.  6)  expressed  in 
terms  of  the  stagnation  enthalpy, 

^  11  ri  11  V*  ri  11^ 

*  =  r  E  +  [1  -  -]  y  +  -  -]  g  hjn,,  (10) 


where  V  is  the  magnitude  of  the  turbulent  (fluctuating)  velocity,  y/k.  The  Schmidt  number, 
Sc  =  peffKpF),  (ratio  of  momentum  transport  to  mass  transport)  is  assumed  to  be  unity. 
The  Prandtl  number,  Pr  =  CpPatfK,  (ratio  of  momentum  transport  to  heat  transport)  is 
assumed  to  be  nearly  unity  (.9)  which  is  considered  adequate  for  gaseous  flows  (Bradshaw 
1981).  The  thermal  conductivity  of  the  gas  mixture  is  denoted  k. 
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The  mixture  equation  of  state  for  a  thermally  perfect  gas  follows  from  Dalton’s  Law, 

(11) 

where  $  =  -Mj  is  the  molecular  weight  of  species  j,  and  X  is  the  specific  gas 

constant.  Equation  11  is  used  to  recover  the  density  from  the  numerical  solution  of  the 
governing  equations. 


The  stream  function-vorticity  form  of  the  governing  equations  has  been  widely  utilized 
2md  faicilitates  the  use  of  numerically  efficient  Gauss-Seidel  relaxation  algorithms.  Stream 
function,  0,  and  vorticity,  w  are  defined  using  (Schlichting  1979), 


The  governing  equations,  derived  by  Nusca  (1989),  can  be  expressed  in  the  form  of  a  general 
variable,  0.  This  variable  can  represent  stream-function,  vorticity,  azimuthal  velocity,  stag¬ 
nation  enthsdpy,  species  mass  fraction,  turbulence  kinetic  energy,  or  turbulence  dissipation 
rate. 


\-l 

dz  \ 

dr  \ 

‘1)1 

dz 

_  A 

dr 

+  rd^  =  0  (13) 


For  example,  0  =  0,a^  =  0,6^  =  l/(pr*),c^  =  l,d^  =  -u//r  (which  yields  Eq.  12).  For 
N  species  only  N  —  I  specie  equations  (0  =  mj)  must  be  solved,  since  the  sum  of  the  mass 
fractions  must  equal  unity.  In  effect  the  global  continuity  equation  (Eq.  1)  is  the  A^th  specie 
equation  since  the  summation  of  all  specie  equations  yields  the  continuity  equation. 


The  pressure  can  be  recovered  from  the  0  —  a;  form  of  the  equations  after  a  converged 
solution  of  Equation  13  is  achieved  or  after  each  iteration,  if  pressure  variations  are  expected 
to  have  a  significant  effect  on  density.  The  radial  and  axial  momentum  equations  (Eq.  2  and 
3)  can  be  rearranged  to  yield: 

~-Pi{r,z)  and  |^  =  -P2(r,z)  (14) 

where  Pi  and  P2  are  functions  of  p,  V  and  f.  Along  any  path  from  point  A  to  point  B  in 
the  flowfield,  the  pressure  is  given  by; 

PB-PA=  I  {P\dz  -f  Pidr)  (15) 

J  A 

Since  p  is  a  scalar,  ps  —  Pa  should  be  path  independent  and  therefore,  can  serve  as  a 
consistency  check  on  the  converged  solution.  In  most  cases  a  pressure  difference  is  desired 
to  form  the  pressure  coefficient  at  a  point.  However,  if  the  pressure  at  point  B  is  required 
the  pressure  at  point  A  is  assigned  to  a  known  inlet  value  and  integration  proceeds  from  the 
inlet  to  point  B  in  the  flowfield. 
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2.2  Turbulence  Model.  A  two-equation  turbulence  model  has  been  suggested  by 
Kim  and  Chung  (1989)  for  multiple  species  flows.  This  model  describes  the  turbulence 
viscosity  (/it)  as  a  function  of  turbulence  kinetic  energy  (k)  and  dissipation  rate  (e)  as 
=  pCsk^fe.  A  set  of  partial  differential  equations  is  written  for  k  and  c. 


dk 


dk  1 


■£ 

( 

d 

(  3t\i 

dz 

k 

4-  - 

^  dr 

de 

% 


de 

% 


1  ra  / 

dey 

CiGe 

;  [dz 

dz ) 

r'arjj 

k 

^dz)  \ 

dr  j 

»  ^[r 

)  j  var 

dz  j 

p  +  Pt/K, 

=  hK 

=  1.3,  C,  = 

1.44, 

C-ipe^ 


(16) 

(17) 

(18) 
.09. 


These  equations  are  solved  along  with  the  Navier-Stokes  equations  (Eq.  13)  using, 

<k  =  k,a^  =  1,6^  =  =  l,d^  =  r{G  -  pe) 

^  =  c,a^  =  1,6^  =  r/i<,c^  =  l,d^  =  r{CiGtlk  -  C^pe^/k) 


2.3  Boundary  and  Initial  Conditions.  The  boundaries  of  the  inner  grenade  cham¬ 
ber  (see  Figure  1)  are  the  inlet  (connected  to  the  pyrotechnic  chamber),  the  exit  (nozzle 
throat),  the  symmetry  axis,  the  chamber  wail  lined  with  payload  material,  and  a  section  of 
solid  wall  along  the  nozzle.  Since  the  governing  equations  (Eq.  13)  are  elliptic  the  conditions 
along  these  boundaries  must  prescribe  values  of  the  dependent  vauriables,  the  gradient  of  the 
dependent  variables  in  the  norm2d  direction,  or  an  algebraic  relation  which  connects  the 
values  of  the  dependent  variables  to  the  normal  component  of  velocity. 

At  the  inlet  plane,  radial  profiles  of  all  dependent  variables,  xl;,uj,mj,k,v,k,  and  c  as 
well  as  values  for  K,  T,  p,  p,  and  are  specified.  It  is  assumed  that  the  flow  at  the 
inlet  plane  consists  of  air  and  that  the  diffusion  of  payload  into  the  airstre^un  from  the 
chamber  wall  does  not  effect  the  inlet  flow.  During  actual  grenade  operation  the  inlet  flow 
consists  of  pyrotechnic  combustion  gases  (e.g.  H3O,  CO3)  the  exact  nature  of  which  was 
unavsulable  for  the  numerical  simulations.  A  subsonic  inlet  flow  velocity  assumption  is  used 
in  accordance  with  the  elliptic  nature  of  the  governing  equations.  Initial  conditions  for  all 
dependent  variables  are  supplied  by  the  inlet  boundary  conditions.  The  turbulence  model  is 
initialized  using  k  =  koo  =  o‘V^,€  =  ®  The  initial  turbulence  kinetic 

energy  is  specified  as  (1/q)%  of  the  inlet  kinetic  energy. 

The  exit  plane  is  located  at  the  nozzle  throat  where  the  flow  is  assumed  to  be  subsonic. 
The  streamlines  at  the  exit  plane  are  assumed  to  be  parallel  to  the  symmetry  axis;  thus,  the 
gradients  of  ail  dependent  variables  along  these  streamlines  are  zero.  These  assumptions  are 
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reasonable  since  experience  for  large  Reynolds  numbers  has  shown  that  the  exact  nature  of 
the  exit  plane  boundary  conditions  has  little  effect  on  the  flowfield  solution  when  convection 
is  significant  (i.e.  large  inlet  mass  flow)  (Bradshaw  et.ad.  1981,  Khalil  et.al.  1975). 

For  mass  continuity,  the  symmetry  axis  is  considered  to  be  a  streamline  of  the  flow,  thus 
0  =  constant.  Along  the  symmetry  axis  r  =  0,  thus  d^fdz  =  dtl}fdr  =  0  via  Equation 
12.  The  value  of  0  along  this\^undary  can  be  determined  from  values  of  adjacent  to 
the  boundary  using  a  one-si(^ed  finite-difference  for  dipfdr  at  the  axis.  From  Equation  12, 
the  boundary  value  for  u?  is  zero.  The  axis  boundary  values  for  the  remaining  dependent 
variables,  are  determined  from  =  0. 


The  no-slip  condition  (u  =  «;  =  0,  and  v/r  =  fi,  where  fl  is  the  wall  spinrate)  is 
applied  to  the  solid  walls.  Therefore  =  constant,  via  Equation  12.  For  convenience  ip  =  0 
is  chosen.  One-sided  finite-differences  for  dtp  {dr  and  dtp/dz  are  used  in  Equation  12  to 
determine  the  wall  value  for  tj.  For  an  inert  wall,  the  normal  gradient  of  all  mass  fractions, 
dmjidn.,  are  set  to  zero.  For  a  wall  with  payload  material,  the  boundary  condition  is  based 
on  the  assumption  that  the  payload  material  is  continually  vaporizing  =  l,m,ir  = 

0,  mntixture  =  0).  The  wall  temperature  is  set  to  the  vaporization  temperature  of  the  payload, 
Tw«ii  =  Tvap.  Figure  2  shows  the  results  of  a  thermogravimetric  analysis  of  yellow  dye  (i.e. 
payload  material)  (Turetsky  1991).  The  percent  weight  loss  of  the  sample  is  plotted  as  a 
function  of  temperature.  Note  that  vaporization  of  the  material  (i.e.  significant  weight  loss) 
occurs  over  a  narrow  temperature  range,  supporting  the  use  of  such  a  boundary  condition. 
The  rate  of  burning  (regression  rate)  on  the  payload  surface  varies  as  a  function  of  position 
along  the  surface  and  is  computed  from  the  temperature  gradient  normal  to  the  wall. 


r  = 


—Kp  dT 


Pp^v*p 


dr 


(19) 


where  «p  and  pp  are  the  thermal  conductivity  and  density  of  the  solid  payload,  and  hy^ 
is  the  heat  of  vaporization  of  a  unit  mass  of  payload.  Values  of  the  thermal  conductivity, 
density,  and  heat  of  vaporization  can  be  determined  for  most  grenade  payload  materials. 


The  boundary  conditions  for  the  payload  surface  are  based  on  the  assumption  of  single 
diffusion,  i.e.  diffusion  of  gaseous  payload  molecules  into  the  airstream  without  diffusion  of 
air  molecules  into  the  payload  material.  Single  diffusion  has  been  studied  by  R.D.  Present 
(1958).  The  general  equation  of  mutual /thermal  diffusion  is  given  by, 


njUj 


nF 


IfM 

dz\n) 


n^Tr  dT 
~T^1^ 


O'  =  1.2) 


(20) 


where  =  molecular  density  of  species  j,  n  —  =  convection  velocity  of  species  j, 

z  =  diffusion  direction,  F  =  mass  diffusion  coefficient  (VA/3,  A  =  the  molecular  mean  free 
path),  and  Fr  =  thermal  diffusion  coefficient  («/(phv«p))-  The  equation  of  “single  diffusion” 
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can  be  obtained  by  assuming  that  air  molecules  (i=l)  are  moving  in  a  boundary  layer  so 
that  uj  ~  0  and  payload  molecules  (j=2)  are  at  rest,  U3  =  0.  In  addition,  the  payload 
molecules  are  closely  psicked  so  that  ni  n2.  As  a  result, 

G^  =  Gii^0,  Gpjiyioad  =  (?j  -(ni  +  ~7^dz 

The  mass  diffusion  terms  of  Equation  21  are  incorporated  in  the  Jj  terms  of  Equation  5  and 
Fx  is  customarily  neglected. 

The  payload  surface  boundary  condition  may  also  be  prescribed  using  a  surface  ablation 
model  such  as  described  by  Moss  (1976).  In  this  model  the  paylo2td  vaporization  temper¬ 
ature  and  heat  of  vaporization  aure  prescribed  as  functions  of  the  pressure  applied  to  the 
material  surface.  An  equation  similar  to  Ekjuation  19  is  used  along  with  these  functions  in  a 
iterative/coupling  procedure  with  the  governing  equations  in  order  to  prescribe  the  payload 
surface  boundary  condition.  The  use  of  an  ablation  model,  as  opposed  to  assuming  that  the 
payload  surface  is  continually  vaporizing  (TVau  =  TUp)  is  more  critical  for  cases  where  the 
vaporization  characteristics  of  the  material  are  a  strong  function  of  pressure  and  temperature 
(i.e.  unlike  that  shown  in  Fig.  2).  The  major  drawback  to  using  an  ablation  model  is  the 
incre<ised  computational  cost  of  the  iterative  procedure  and  the  requirement  of  experimental 
data  in  the  form  7vap  ™  ^vap(p)7  hvap  —  ^vap(p)‘ 

Boundary  conditions  for  turbulence  variables,  k  and  e,  are  A:  =  0,  c  =  .056/i(5u/5r)*//> 
for  solid  walls  and  k  =  10“®V^,£  =  for  the  inlet  flow.  Along  the  chamber  axis  and 

exit  plane  dkfdr  =  de/dr  =  O^dk/dz  =  dt/dz  —  0,  respectively. 


2.4  Computational  Algorithm.  Equation  13  can  be  reduced  to  a  successive-substitution 
formula  for  flow  variable  <f>  at  each  node  on  the  computational  grid.  Central  finite-differences 
are  used  for  the  diffusive  and  source  terms  and  upwind  differences  for  the  convective  terms. 
Using  upwind  differencing  in  the  species  conservation  equations  (E)q.  5)  reduces  the  occur¬ 
rence  of  negative  species  mass  fractions  in  mixing  layers.  The  resulting  system  of  equations 
for  the  entire  grid  is  solved  using  a  Gauss-Seidel  relaxation  scheme  (Nusca  1989).  Each  iter¬ 
ation  cycle  is  made  up  of  M  sub-cycles,  where  M  is  the  number  of  equations  being  considered 
(M  must  be  at  least  2  since  the  equations  for  <f>  =  u;fr  and  ^  =  0  2ure  the  minimum  required 
to  define  the  flow).  In  each  sub-cycle,  grid  points  are  scanned  row  by  row  and  a  single 
variable  is  updated.  The  variables  (j/r  and  ^  are  updated  in  order  followed  by  all  other 
variables.  When  all  sub-cycles  are  completed  a  new  iteration  cycle  is  started  in  which  the 
values  of  the  variables  from  the  latest  iteration  are  immediately  used.  This  is  consistent  with 
the  Gauss-Seidel  methodology  (Carnahan  1969).  Convergence  is  satisfied  when  the  greatest 
relative  change  in  any  flow  variable  is  smaller  than  a  prescribed  tolerance. 
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3.  RESULTS 


In  order  to  demonstrate  the  numerical  method,  simulations  were  performed  for  an  exper¬ 
imental  grenade  model  loaded  with  yellow  dye  payload  simulant  that  has  been  the  subject 
of  testing  at  CRDEC.  In  these  tests  the  grenade’s  inner  chamber  was  instrumented  with 
pressure  transducers  and  thermocouples  (for  temperature  measurements).  The  chamber  was 
lined  with  yellow  dye  payload  material  about  .25  inches  in  thickness.  A  typical  test  run  yields 
measured  values  for  the  chamber  inlet  pressure  and  temperature  (due  mainly  to  pyrotechnic 
ignition)  of  about  3.5  psig  (1.24  atm)  and  425  C.  The  nozzle  or  exit  plane  values  were  mea¬ 
sured  as  1.5  psig  (1.1  atm)  and  375  C.  These  values  were  used  as  boundary  conditions  for 
the  numerical  simulation.  The  chemical  properties  of  the  yellow  dye  payload  material  were 
measured  as  Tv^p  =  241  C,  Av*p  =  102.9  J/g,  A4  =  273  g/mole.  The  density  and  thermjd 
conductivity  of  the  dye  have  not  been  measured  but  were  taken  as  pp  =  1.8  g/mole  and  Kp  = 
.00143  cal/s-cm-C,  respectively.  Sensitivity  studies  using  the  present  numerical  model  show 
that  the  predicted  payload  regression  rate  depends  significantly  on  the  value  of  /Cp  and  less 
on  the  value  of  Pp.  The  inlet  gas  was  assumed  to  be  air  with  M  —  28.8  g/mole  and  the  M 
of  the  air/payload  gas  mixture  taken  as  150  g/mole. 

Figure  3  shows  the  computational  grid  (or  mesh)  used  for  flowfield  simulation  within 
the  grenade’s  inner  chamber,  bounded  axially  by  the  inlet  and  nozzle  and  radially  by  the 
symmetry  axis  and  chamber  wall.  The  chamber  is  cylindrical  with  an  overall  length  of  4.0 
inches  aind  .375  inches  in  diameter.  These  dimensions  represent  the  chamber  shape  before 
payload  surface  regression.  Note  that  the  vertical  axis  in  Figure  3  has  been  magnified  for 
clarity  by  a  factor  of  10  over  the  horizontal  axis.  About  131  grid  points  were  used  in  the  axial 
direction  with  35  points  in  the  radial  direction;  the  grid  points  were  unevenly  distributed  in 
order  to  cluster  points  along  the  boundaries  (except  at  the  symmetry  axis).  Figure  4  shows 
the  computed  velocity  vectors  (i.e.  arrows  whose  length  is  representative  of  the  magnitude 
of  local  gas  velocity  and  direction  indicates  the  orientation  of  the  local  velocity)  displayed  at 
every  6th  axial  grid  point  and  every  other  radial  point  (for  clarity).  A  thick  boundary  layer 
that  develops  along  the  chamber  wall  (i.e.  payload  surface)  can  be  observed.  The  computed 
average  chamber  exit  velocity  is  approximately  121  ft/s.  The  surface  of  the  payload  material 
was  allowed  to  ablate  (i.e.  regress)  for  .25  seconds  resulting  in  the  shape  shown  in  Figure  5 
(computational  grid).  Note  that  the  material  surface  has  been  blunted  at  the  chamber  inlet 
but  more  evenly  eroded  over  most  of  the  surface  with  the  exception  of  a  discontinuity  at  the 
payload/nozzle  wall  junction  (z  =  3.5  inches).  Figure  6  shows  the  computed  velocity  vectors 
for  this  chamber  shape.  Due  to  the  contour  of  the  chamber  entrance,  the  flow  is  significantly 
accelerated  (indicated  by  longer  arrows  in  the  figure)  forming  a  thinner  boundary  layer 
along  the  payload  surface  (i.e.  Izu'ger  Re).  The  computed  average  chamber  exit  velocity  is 
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approximately  276  ft/s.  Figure  7  shows  the  distribution  of  payload  regression  rate  (Eq.  19) 
along  the  payload  surface  for  both  the  initial  time  and  .25  seconds  elapsed  time.  Payload 
regression  rate  is  largest  near  the  chamber  inlet  and  nearly  uniform  over  most  of  the  payload 
surface.  The  regression  rate  increases  with  elapsed  time. 

Figure  8  shows  radial  profiles  of  2ixial  velocity  component  at  both  the  initial  and  elapsed 
time  as  well  as  two  axial  positions  along  the  chamber  -  midlength,  z  —  2in,  and  the  exit 
plane,  z  =  4in.  Consistent  with  the  imposed  boundary  conditions  the  velocity  is  zero  at  the 
chamber  wall  (r  =  .1875in)  and  the  velocity  gradient  is  zero  at  the  chamber  axis  (r  =  0). 
The  retarding  effect  of  the  wall  boundary  layer  on  the  velocity  profile  can  also  be  observed. 
Note  that  the  profiles  at  chamber  midlength  show  the  payload  surface  regressing  from  r  = 
.1875in  to  about  .195in  from  the  centerline  while  the  profiles  at  chamber  exit  show  no  wall 
regression  since  the  chamber  wall  is  solid  at  this  location.  Comparing  the  chamber  midlength 
profiles  with  those  at  the  exit  plane  show  that  the  boundary  layer  (or  mixing  layer)  thickens 
with  axial  location  downstream  of  the  inlet.  Comparing  the  initial  time  profiles  with  those 
at  the  elapsed  time  show  that  the  gas  flow  is  significantly  accelerated  as  the  payload  surface 
regresses  and  the  boundary  layer  thickness  decreases  (i.e.  larger  Re). 

Figure  9  shows  radial  profiles  of  gas  temperature  at  both  the  initial  and  elapsed  time 
as  well  as  two  axial  positions  along  the  chamber  -  midlength,  z  =  2in,  and  the  exit  plane, 
z  =  4in.  Consistent  with  the  imposed  boundary  conditions  the  temperature  is  Ty^,  at  the 
payload  surface  while  the  temperature  gradient  is  zero  at  both  the  chamber  axis  and  the 
solid  chamber  wall  at  the  nozzle  (adiabatic  wall  condition).  Note  that  the  profiles  at  chamber 
midlength  show  the  payload  surface  regressinj  from  r  —  .1875in  to  about  .195in  from  the 
centerline  while  the  profiles  at  chamber  exit  show  no  wall  regression  since  the  chamber  wall 
is  solid  at  this  location.  Within  the  boundary  layer  established  on  the  payload  surface  the 
temperature  of  the  air/payload  mixture  is  gradually  increased  until  it  reaches  the  centerline 
(i.e.  inlet)  value.  Comparing  the  chamber  midlength  profiles  with  those  at  the  exit  plane 
show  that  the  thermal  boundary  layer  thickens  with  axial  location  downstream  of  the  inlet. 
Comparing  the  initial  time  profiles  with  those  at  the  elapsed  time  show  that  the  gas  flow  is 
slightly  hotter  at  any  chosen  radial  position  within  the  mixing  layer. 

Figures  10  and  11  show  radial  profiles  of  mixture  (payload  +  air)  mass  fraction  at  both 
the  initial  and  elapsed  time  as  well  as  two  axial  positions  along  the  chamber  -  midlength, 
z  =  2in,  and  the  exit  plane,  z  =  4in.  Consistent  with  the  imposed  boundary  conditions 
the  mixture  mass  fraction,  mnu*ture>  is  zero  at  the  payload  surface  (i.e.  mp*yio»d  =  1, 

=  0)  and  along  the  chamber  centerline  (i.e.  mp^yioad  =  0,  mtir  =1).  Along  solid  walls  (i.e. 
chamber  exit)  and  at  the  chamber  centerline  the  gradient  of  m„jxture  is  zero.  The  mixture 
mass  fraction  is  exactly  or  nearly  unity  at  some  point  in  the  flowfield  where  mpayioAd  and 
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m*ir  are  equivalent.  Note  that  the  profiles  at  chamber  midlength  show  the  payload  surface 
regressing  from  r  =  .1875in  to  about  .195in  from  the  centerline  while  the  profiles  at  chamber 
exit  show  no  wall  regression  since  the  chamber  wall  is  solid  at  this  location.  Comparing 
the  chamber  midlength  profiles  with  those  at  the  exit  plane  show  that  the  chamber  core 
flow,  consisting  of  air  (mmixture  =  0),  thins  from  a  radial  position  of  about  .I2in  to  .lOin  at 
the  initial  time  and  .15in  to  .13in  at  the  elapsed  time.  The  diffusion  of  payload  material 
into  the  core  flow  is  advanced  with  axial  position  as  the  boundary  layer  imposed  on  the 
payload  surface  grows.  Comparing  the  initial  time  profiles  with  those  at  the  elapsed  time 
show  that  the  mixing  layer  (i.e.  highest  mixture  concentration)  follows  the  payload  surface 
as  regression  progresses. 

4.  CONCLUSIONS 

Material  dissemination  from  the  payload  chamber  inside  a  dual-chamber  grenade  has  been 
simulated  using  computational  fluid  dynamics.  A  thermal  dissemination  technique  that  uses 
the  hot,  moving  gases  generated  from  combustion  of  a  pyrotechnic  within  the  grenade  has 
been  investigated.  The  dissemination  process  is  initiated  by  ablation  and  vaporization  of 
the  payload  material  from  exposure  to  the  thermal  effects  of  pyrotechnic  combustion.  This 
material  in  gaseous  form  is  entrained  in  this  flow  and  expelled  from  the  grenade.  The  Navier- 
Stokes  equations  along  with  chemical  species  conservation  equations  were  used  to  simulate 
the  diffusion  and  convection  processes  of  the  flowfield  within  the  grenade’s  payload  chamber. 
Numerical  simulations  reveal  that  diffusion  of  the  payload  material  is  accomplished  within 
a  boundary  layer  that  is  established  along  the  payload  surf£u:e  (chamber  wall)  and  that  a 
core  flow,  basically  unaffected  by  the  payload  ablation,  resides  over  about  one-half  of  the 
chamber  diameter.  Thermal  exposure  of  the  payload  mixture  to  the  hot  combustion  gases 
from  the  pyrotechnic  is  concentrated  in  this  layer  where  the  temperature  is  below  that  of  the 
core  flow  gases.  As  the  payload  surface  regresses,  the  flow  thru  the  chamber  is  accelerated 
and  the  diffusion/mixing  layer  follows  the  regressing  surface  which  results  in  an  expansion 
of  the  core  flow.  The  regression  rate  increases  with  elapsed  time  (since  pyrotechnic  ignition) 
as  the  shape  of  the  payload  surface  (chamber  wall)  is  contoured  by  material  ablation. 
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FUNCTIONING 


Figure  1.  Schematic  showing  operation  of  the  dual-chamber  grenade  thermal 

dissemination  (M.  Miller,  U.S.  Army  CRDEC,  used  with  permission) 
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Figure  5.  Computational  grid  for  cheunber  shape  after  .25  seconds  of  payload  surface 
ablation 
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Figure  8.  Computed  radial  profiles  of  axial  component  of  gas  flow  velocity;  results  for 
initial  and  elapsed  time  shown  &t  the  chamber  midlength  and  nozzle  exit 
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Figure  9.  Computed  radial  profiles  of  gas  temperature;  results  for  initial  and  elapsed  time 
shown  at  the  chamber  midlength  and  nozzle  exit 
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Figure  10.  Computed  radial  profiles  of  gas  composition  (mixture  mass  fraction);  results 
for  initial  and  elapsed  time  shown  at  the  chamber  midlengtb 


Figure  11.  Computed  radial  profiles  of  gas  composition  (mixture  mass  fraction);  results 
for  initial  and  elapsed  time  shown  at  the  chamber  nozzle  exit 
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LIST  OF  SYMBOLS 


Cp  specific  heat  capacity,  constant  p 

h  molar  specific  enthalpy 

k  total  enthalpy 

fivap  payload  material  heat  of  vaporization 
J  flow  rate  or  flux 

k  turbulence  kinetic  energy 

L  characteristic  length 

m  species  mass  fraction 

M  moleculau'  weight 

n  molecular  density 

N  number  of  species 

p  static  pressure 

Pr  Prandtl  Number 

r  radial  direction 

r  payload  surface  regression  rate 

31?  specific  gas  constant,  (7  —  l)cp/7 

R  universal  gas  constant,  3?  Mj 

Re  Reynolds  Number 

Sc  Schmidt  Number 

t  time 

T  static  temperature 

u  radial  velocity  component 

V  azimuthal  velocity  component 

V  magnitude  of  the  local  velocity  vector 

V  uf  +  v0  +  wz 

w  axial  velocity  component 

X  mole  fraction 

z  axial  direction 
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Greek  Symbols 


7  ratio  of  specific  heats 

r  diffusion  coefficient 

(  turbulence  dissipation  rate 

$  azimuthal  direction 

K  thermal  conductivity 

fi  molecular  viscosity 

p  density 

T  shear  stress  vector 

^  general  flow  variable 

^  stream  function 

u;  vorticity 

Superscripts 

unit  vector 

total  or  stagnation  quantity 
rate 


Subscripts 

eff  effective 

h  enthalpy  contribution 

j  j-th  mixture  component  or  species 

k  turbulence  kinetic  energy  contribution 

p  constant  pressure 

p  payload  quantity 

r  radial  component  or  radial  direction 

t  turbulence  quantity 
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z  axial  component 

c  turbulence  dissipation  rate 

$  azimuthal  component 

oo  freestream  quantity 
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East  Amherst,  NY  14501-0305 

1  Universal  Propulsion  Company 

ATTN:  H.J.  McSpadden 
25401  North  Central  Ave. 

Phoenix,  AZ  85027-7837 

1  SRI  International 

Propulsion  Sciences  Division 
ATTN:  Tech  Library 
333  Ravenwood  Avenue 
Menlo  Park,  CA  94025-3493 


No.  of 

Copies  Organization 


1  Cdr,  USACSTA 

ATTN:  STECS  PO/R.  Hendricksen 
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No.  of 

Coates  Organization 

1  Ernst-Mach-lnstitut 

ATTN:  Dr.  R.  Reiser 
Haupstrasse  18 
Weil  am  Rheim 
Germany 

1  Defence  Research  Agency.  Military 

Division 

ATTN:  C.  Woodley 
RARDE  Fort  Halstead 
Sevenoaks.  Kent,  TNI  4  7BP 
England 

1  School  of  Mechanical,  Materials,  and 

Civil  Engineering 
ATTN;  Dr.  Bryan  Lawton 
Royal  Military  College  of  Science 
Shrivenham.  Swindon,  Wiltshire. 

SN6  8LA 
England 


No.  of 

Copies  Organization 

1  Explosive  Ordnance  Division 

ATTN;  A.  Wiidegger-Gaissmaier 
Defence  Science  and  Technology 
Organisation 
P.O.  Box  1750 

Salisbury.  South  Australia  5108 

1  Armaments  Division 

ATTN:  Dr.  J.  Lavigne 
Defence  Research  Establishment 
Valcartier 

2459.  Pie  XI  Blvd..  North 
P.O.  Box  8800 

Courcelette,  Quebec  GOA  1 RO 
Canada 

1  U.S.  Army  European  Research  Office 

ATTN;  Dr.  Roy  E.  Richenbach 
Box  65 

FPO  New  York  09510-1500 


2  Institut  Saint  Louis 

ATTN:  Dr.  Marc  Giraud 

Dr.  Gunther  Sheets 
Postfach  1260 
7858  Weail  am  Rhein  1 
Germany 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


TOs  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your 
comments/answers  to  the  iicms/questions  below  will  aid  us  in  our  efforts. 

1 .  ARL  Report  Number  ARL-TR-77 _ Date  of  Report  February  1993 _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for 

which  the  report  wiU  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of 
ideas,  etc.) _ _ _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved, 
operating  costs  avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate 
changes  to  organization,  technical  content,  format,  etc.) _ 


Organization 


CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address 
above  and  the  Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  staple  or  tape  closed,  arxl  mail.) 
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OFRCIAL  BUSINESS 


BUSINESS  REPLY  i\L4IL 

RRST  Ct>SS  URMT  No  00)1,  APS.  MD 

Postage  *ill  Oe  paid  by  addressee 


Director 

U.S.  Army  Research  Laboratory 
ATTN:  AMSRL-OP-CI-B  (Tech  Lib) 
Aberdeen  Proving  Ground,  MD  21005-5066 


